Neutrino mass in cosmology: status and prospects 



Yvonne Y. Y. Wong 

Institut fiir Theoretische Teilchenphysik und Komsologie, RWTH Aachen 

D-52056 Aachen, Germany 
yvonne.wong@physik.rwth-aachen.de 

November 8, 2011 
Abstract 

I give an overview of the effects of neutrino masses in cosmology, focussing on the 
role they play in the evolution of cosmological perturbations. I discuss how recent ob- 
servations of the cosmic microwave background anisotropies and the large-scale matter 
distribution can probe neutrino masses with greater precision than current laboratory 
experiments. I describe several new techniques that will be used to probe cosmology 
in the fiiture, as well as recent advances in the computation of the nonlinear matter 
power spectrum and related observables. 
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1 Introduction 



Standard big bang theory predicts some 10^^ neutrinos per flavour in the visible universe, an 
abundance second only to the cosmic microwave background (CMB) photons. Relativistic 
neutrinos constitute a significant fraction of the total energy density in the epoch of radiation 
domination, and play an important role in the primordial synthesis of light elements. Even 
in today's matter and dark energy-dominated universe, a minimum neutrino mass of 0.07 eV 
as established by the solar and atmospheric neutrino oscillations experiments [1] means that 
relic big bang neutrinos can still contribute at least 0.5% of the total matter density. These 
neutrinos have a definite impact on cosmic structure formation. 

This article provides a brief overview of the subject of neutrino cosmology, focussing on 
how precision cosmological probes such as CMB anisotropics and the large-scale structure 
(LSS) distribution of matter can be used to determine or constrain the absolute neutrino 
mass scale. Readers interested in the role of neutrinos for big bang nucleosynthesis should 
consult the review of [2]. After introducing the relic neutrino background in section 2, I 
outline the theoretical framework for computing precision cosmological observables using 
linear perturbation theory in section 3. In section 4 I describe the capacities of present 
observations to probe the absolute neutrino mass scale, while in sections 5 and 6 I report 
on some recent developments in the computation of nonlinear corrections to the LSS matter 
power spectrum and the halo mass function. Section 7 contains descriptions of various future 
observations and their forecasted sensitivities to the neutrino mass. I conclude in section 8. 

2 The relic neutrino background 

2.1 Thermal history of the neutrino: decoupling 

In the first second after the big bang, the temperature of the universe is so high that even 
weak interactions are able to hold neutrinos in thermal equilibrium with the cosmic plasma 
(photons, electrons, positrons, etc.). When in equilibrium, the neutrinos' phase space density 
follows the Fermi-Dirac distribution. 



where the neutrinos share a common temperature with the cosmic plasma Ty = = T, 
and the chemical potential ji is expected, on theoretical grounds, to be similar to the uni- 
versal matter-antimatter asymmetry, /^/Tj, ~ 10"^'^, i.e., negligible. The number and energy 
densities per flavour are given respectively by 



where g = 2 denotes the two helicity states, i.e., the expressions (2) count both neutrinos 
and antineutrinos. 




(1) 
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As the universe expands and cools, the weak interaction rate drops off with temperature 
as r ~ G\T^. When interactions become too infrequent to match the Hubble expansion 
H ~ T^/mpiancki neutrinos decouple from the cosmic plasma. By demanding F ~ i7, we see 
that neutrinos decouple at T ~ 1 MeV. Laboratory experiments indicate that the absolute 
neutrino mass should not exceed 0{1) eV [3,4]. Thus, at decoupling, the neutrinos are 
relativistic, and their phase space density is well- approximated by the relativistic Fermi- 
Dirac distribution, 

UE,T,)^ \ (3) 

1 + exp[p/T^J 

where for simplicity the chemical potential has been set to zero. 



2.2 After decoupling: the neutrino phase space and temperature 

Because neutrinos decouple while they are relativistic, they preserve, to an excellent approx- 
imation, their relativistic Fermi-Dirac phase space density even after they become nonrel- 
ativistic at later times. The neutrino phase space after decoupling is a frequent point of 
confusion and deserves a detailed explanation. 

Consider neutrinos with mass m,j = 1 eV decoupling at a temperature of 1 MeV. By 
integrating the Fermi-Dirac distribution (1), one finds that only 1 in 10^^ neutrinos has a 
momentum p < rUi,. This justifies the use of equation (3) to describe the neutrino phase 
space distribution up to the point of decoupling. 

Physical momentum p redshifts as a"^, where a is the scale factor, for both massless 
and massive particles. It is thus convenient to define a non-redshifting momentum variable 
y = ap, so that the neutrino number density (2) can be written as oc J dy^ f,y{E,T,j). 
After decoupling, particle number conservation requires that the riu scales as riu oc a~^. 
Thus comparing these two expressions one sees immediately that agreement can only be 
achieved if the phase space density f^{E,T^) is independent of a at all times. Inspection 
of equation (3) shows that this can be brought about if, similar to the physical momentum 
p, the neutrino temperature T„ also redshifts as a~^. Thus the relativistic Fermi-Dirac 
distribution (3) remains an apt description of the neutrino phase space density, even after 
the neutrinos become nonrelativistic. 

Immediately after decoupling the neutrinos have the same temperature as the cosmic 
plasma. This correspondence remains true for some time because both T^y and scale 
as a~^. At T ~ me/3 ~ 0.2 MeV, electrons and positrons become nonrelativistic. Their 
annihilation into photons leads to the transfer of entropy originally residing in the e+e~ fluid 
to the photon fluid, which causes the latter to heat up. This reheating is, however, not felt 
by the decoupled neutrinos. Hence the neutrinos will emerge from the e'^e' annihilation 
event colder than the photons. One can show using entropy conservation arguments (see, 
e.g., [5]) that the new neutrino temperature is related to the photon temperature via 

T.= y T,- (4) 
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This expression remains valid today. 

Equation (4) provides a useful estimate for the temperature of the relic neutrino back- 
ground. However, we do expect some small corrections owing to the fact that the processes 
of neutrino decoupling and e^e^ annihilation occur in close proximity in time and that these 
processes are not instantaneous; At the time of e+e~ annihilation, some neutrinos, particu- 
larly those at the high energy tail of the Fermi-Dirac distribution, are still coupled to the 
cosmic plasma and will partake in the reheating process. Thus in general we expect the 
neutrino energy density to be a little higher than is implied by the relation (4). 

This increase in the neutrino energy density is usually parametcriscd in terms of an in- 
crease in the effective number of neutrino families iVes, defined via X)j p^^i = N^s x p^fi, where 
J2i Pu,i is the total energy density residing in all neutrino species, and pi,^ = {7 / 8) {n"^ /SO) gT^ 
denotes the "standard" neutrino energy density per flavour. Evidently, this relation is 
uniquely defined only at early times when the neutrinos are still rclativistic. Taking also into 
account neutrino flavour oscillations and flnite temperature QED effects, A^eflf = 3.046 [6]. 



2.3 Properties of the relic neutrino background today 

Using equation (4) and measurements of the present CMB temperature Tcmb.o = 2.725 ± 
0.001 K [7] leads us to expect a relic neutrino background of temperature T)^ o = 1.95 K ~ 
10~"^ eV. The number density is expected to be 112 cm~^ per flavour from equation (2). 

The exact energy density per flavour depends on whether the neutrinos are relativistic or 
nonrelativistic today. The neutrinos are relativistic if rrii, <C T„fi, in which case their energy 
density per flavour is p^ = (7/8)(4/ll)'^/^pcMB- Normalised to the present-day critical density 
Pcrit,o = SHq/StiG, where Hq = 100 h is the present Hubble parameter and G Newton's 
constant, we flnd fli,h^ = {pufl/ Pcrit,o)h'^ = 6 x 10^^. In other words, the energy density due 
to relativistic neutrinos today is completely negligible. 

However, if rrii, ^ Ti,^q, then the energy density per flavour is p^ — mi,nv, or, equivalently. 

Thus, even for a neutrino mass as small as rrii, = 0.05 eV, we expect to flnd a non-negligible 
Vli, ~ 0.1%; These neutrinos then form a dark matter component in the universe. By de- 
manding that massive neutrinos not overdose the universe, i.e., Vt^ one can immediately 
set an upper bound on the sum of the neutrino masses, ^niy < 93 eV [8,9]. Historically 
this is flrst upper bound on the neutrino mass from cosmology and is sometimes known as 
the "closure" bound. 



3 Linear cosmological perturbations 

Neutrino dark matter satisfying the closure bound cannot constitute all of the dark matter 
content of the universe because thermal relics that decouple when relativistic come with a 
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large velocity dispersion. Using the relativistic Fermi-Dirac distribution (3), we find a rough 
estimate of this velocity dispersion: 

(^^thermal) - 81(1 + z) ( — ) km 8"^ (6) 

For a ~ 1 eV neutrino, (I'thermai) — 100 km s~^ is comparable to the typical velocity dispersion 
of a galaxy. For dwarf galaxies, the velocity dispersion is even smaller, ~ 0(10) km s~^. Thus 
the relic neutrinos have much too much thermal energy to be squeezed into small volumes 
to form the smaller structures we observe today [10]. In contrast, cold dark matter (CDM) 
has by definition (^thermal) = 0, and is thus not subject to these constraints. 

Nonetheless, even if relic neutrinos cannot form the bulk of the cosmic dark matter, 
because their kinematic properties are so different from those of CDM, their presence at 
even the fi^ ~ 0.1% level must leave a signature in the large-scale cosmological observables. 
Detecting this signature will then allow us to establish the absolute neutrino mass scale via 
equation (5). 

In this section, I outhne the theoretical framework for predicting the effects of massive 
neutrinos on the CMB anisotropics and LSS matter power spectrum via linear perturbation 
theory. For more detailed discussions of linear cosmological perturbation theory in general, 
see, e.g., [11,12]. 

3.1 The homogeneous universe 

The observed universe appears to be homogeneous and isotropic on scales of (9(100) Mpc. 
On these scales space also appears to be expanding. The simplest spacetime metric that 
captures these observational features has the form 

ds^ = g^.^dx^dx" = (r) [-cir^ + -f^^dx^dx^] , (7) 

where r is the conformal time, and = x are the comoving coordinates. The metric (7), 
known as the Friedmann-Lemaitre-Robertson-Walker (FLRW) metric, forms the basis of 
modern cosmology. 

The spatial part of the FLRW metric encodes the local geometry of space, which can 
be (i) fiat and Euclidean, (ii) spherical (i.e., with positive curvature), or (iii) hyperboloid 
(i.e., with negative curvature). Currently, there is no observational evidence for spatial 
curvature [13]. From a theoretical perspective, it is also difficult to reconcile spatial curvature 
with inflationary cosmology (see, e.g., [14]). We therefore consider only the case of flat spatial 
geometry, so that = 

The energy content of the universe is encoded in the stress-energy tensor T^j,. Homo- 
geneity and isotropy imply that there is only one sensible choice, 

n = n = diag(-p, p,p,p), (8) 

where p and p are the spatially averaged energy density and pressure, respectively, of a 
comoving fluid in its rest frame. Expression (8) can be easily generalised to the multi-fluid 
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case with the replacements p — > Pa and p — > Z^aPa, where a sums over all fluids present: 
this usually means CDM a — neutrinos baryons h (i.e., protons and electrons), and 
photons 7. 

Inserting the FLRW metric (7) and the stress-energy tensor (8) into the Einstein equa- 
tion, R^j^y — \g^pR = SnGTi^u, we find the Friedmann equations, 

n = y-J = P, ^ = ^a(p + 3p), (9) 

with • = d/dr. Likewise, conservation of energy-momentum, V^T'^'^ = 0, imphes ^ + 3^(^-1- 
p) = 0, from which we deduce that p oc for nonrelativistic matter, oc for radiation, 
and so forth. 

3.2 Metric perturbations and the perturbed Einstein equation 

Let us now consider the case of small perturbations around the FLRW spacetime. In its 
most general form the perturbed metric can be written as 

ds^ = a^(r){-{l + 2A)dr^ - 2Bidrdx' + [(1 + 2//i)(5y + 2hij]dx'dx^}. (10) 

Here, the quantities A and represent two scalar degrees of freedom (d.o.f.) under rotation; 
Bi is a vector field which can be further decomposed into a curl-free scalar part (1 
d.o.f.) and a divergence-free vector part B^'"^ (2 d.o.f.); hij is a traceless 3x3 matrix 
decomposable into a scalar /longitudinal h^f (1 d.o.f.), a vector/solenoidal h\f (2 d.o.f.) and 

a tensor/transverse hfj (2 d.o.f) component. Starting with the general expression for the 
stress-energy tensor, 

J^Hi^ ^ fHU ^ ^^^^ ^^^i,^u ^ ^g^v ^ ^^v^ 

where u^^ = dx^ / {—ds^Y^"^ is the fiuid's 4-velocity, p and p the rest frame energy density 
and pressure respectively, and S the shear stress, we can likewise decompose 5T^y under the 
rotation group. The advantage of this decomposition is that once we insert the perturbed 
metric and stress-energy tensor into the Einstein equation, the resulting evolution equations 
for the scalar, vector and tensor perturbed components are completely decoupled from one 
another at linear order. 

Of particular interest to us are the scalar modes. The scalar part of the perturbed 
stress-energy tensor takes the form 

5T\^-5p, 5T\^-{p + p)v^^)\ 

5T\ = {p + p)v^^ - b\'\ ST) = Sp Si + (12) 

where p = p + Sp = p{l + S) and similarly for p, and ff*'' denotes the scalar component of 
the coordinate velocity Vi = dxi/dr. The decoupling of the rotational modes at the linear 
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level means that we are free to consider only the scalar part of the metric perturbations, i.e., 
instead of handling simultaneously 10 d.o.f. of the general perturbed metric (10), we need 
only to deal with the 4 scalar d.o.f.. A, Hl, b\^^ and /i-j^ Of these 4 d.o.f., only two are 
physical; the others are gauge modes related to our choice of the background coordinates. 
In other words, we are free to set any 2 d.o.f. to zero. Doing so is termed fixing or choosing 
a gauge. 

A number of popular gauges are used today in cosmological perturbation theory, each 
with its own advantages (e.g., physical intuition, numerical stability, etc.). I will use the 
conformal Newtonian gauge, which consists of setting i?}*'* = and h!fj — 0, so that the 
scalar perturbed FLRW metric now reads 

ds^ = a^{T){-[l + 2V'(x, T)]dr^ + [1 - 20(x, T)]5ijdx'dx^}, (13) 

where we have also relabelled A = ip and Hl = —0. 

One further decomposition of convenience is the Fourier decomposition, given that in 
linear theory each eigenmode of the Laplacian evolves independently from each other. Define 
the Fourier transform as 

i(k, r) = ^[^(x, r)] = / ^^(x, r)e-^''-. (14) 

We can now proceed to write down the perturbed Einstein equation for the scalar modes in 
Fourier space to linear order in the perturbations: 

4) + ni^ = A7rGa'^(p + p)e, 
4> + H{i^ + 2(t)) + {2H + H^)ij + —{(t)-i^) = inGa'^Sp, 

A;2(0_^) = l27rGa^{p + p)a, (15) 

where /c = |k|, ^ = ik'^Vi is the velocity divergence, and {p + P)(t = — {k~^kikj — |(5jj)E*-' the 
anisotropic stress. Note that these constructions naturally project out the scalar components 
of Vi and E*^-. 

3.3 Boltzmann equation 

Now it remains to write down the evolution equations for the quantities appearing on the 
right hand side of the perturbed Einstein equation (15). The starting point is the phase space 
density of a fluid f{x\Pj,T), defined so that dN = g / {2nY' f{x\ Pj,T)dx^dx'^dx^dPidP2dP^ 
gives the number of particles in a phase space volume dx^dx'^dx^dPidP^dP^. Here, is the 
4-momentum of a fluid element, related to the proper momentum = p measured by a 
comoving observer in the coordinates (13) via Pi = a{l — (f))pi. It is also useful to define a 
new momentum variable q = ap, since in the case = 0, g = |q| is non-redshifting. 
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The general covariant expression for the stress-energy tensor in terms of the phase space 
density of a fluid /(x*, Pj, r) is 

9 f JT-, jn JD / .\-l/2-PM-Pf 

Thus, in terms of the new momentum variable q, we find for the components of the stress- 
energy tensor to first order in the metric perturbations 

(Pq 
(Pq 

T\ = a-V^^/(x,q,r), (17) 



n = -a-^/^6(g,r)/(x,q,r), 



' J (27r)3 . 

where e{q, r) = aE = \fWmF^^~(f' ^ and E is the energy measured by a stationary observer 
with respect to the coordinates (13). 

The phase space density /(x, q, r) can be tracked with the Boltzmann equation, 

where the term C\f \ encapsulates the effects of non-gravitational interactions. Because we are 
interested in the evolution of perturbations in the neutrino fiuid after they decouple from the 
cosmic plasma, this term can be neglected. Equation (18) is nonlinear. Wc therefore resort to 
perturbation theory, and expand the phase space density as /(x, q, r) = /o(q') + /i(x, q, r) + 
■ ■ where /o(q') is the phase space density of the fluid if it were in a homogeneous and 
isotropic FLRW background. For a neutrino fluid, the homogeneous and isotropic component 
of its phase density is /o(5) = [1 + exp(g/T,^o)]~\ where T^^ is the neutrino temperature 
today, (cf. equation (3)). Then, expanding out the Boltzmann equation (18), we find for the 
first order correction, 

^ + S.^ + ^^^0, (19) 
dr e dx. dr dq 

or, equivalently in Fourier space, 

^ + ^^(k . q)A + [qi - ie(k • q)^]^ = 0, (20) 

where q = q/g, and the expression for dq/dr arises from solving for the metric (13) the 
geodesic equation, dP^'/dr = —{l/m)T'^^P/^P'^. 

Consider first the case of CDM. Generally we are more interested in bulk quantities than 
in the detailed phase space distribution. In the case of CDM, two essential bulk quantities 
are the density contrast and the velocity divergence. 
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where n = (27ra)~^/ d^qfo is the mean CDM number density. Evidently, Sc and 9c are 
simply kinetic moments of the phase space distribution. We can therefore construct equa- 
tions of motion for these quantities directly by taking kinetic moments of the Boltzmann 
equation (20) itself. The resulting equations for the zeroth and first moments are 

+ - 30 = 0, Oc + Wc - k'^i^ = 0, (22) 

which are none other than the linearised continuity and Euler equations respectively. Note 
that in order to write the Euler equation in this form, we have had to neglect those phase 
space integrals containing powers of (?/e)" with n> 1. This is a justified approximation for 
CDM since CDM is by definition nonrelativistic. By the same argument, all higher kinetic 
moments are negligible. 

For neutrinos, these approximations are clearly ill-justified, especially at early times 
when they are relativistic. Even at late times it is important that we keep track of kinetic 
moments beyond the first, because it is precisely these higher moments that give rise to 
the non-clustering effects! These considerations inevitably lead us to the conclusion that we 
must solve the full Boltzmann equation (20) for a neutrino fluid. This can be done by brute 
force integration of the equation for a large sample in momentum space. A cleverer route, 
however, is to recast equation (20) into a so-called Boltzmann hierarchy. 

The Fourier space phase space density /i is nominally a function of the wavevector k, 
momentum q, and conformal time r. However, an inspection of equation (20) reveals that 
aside from the dependence on the magnitudes k and fi has no explicit dependence on the 
direction of either vector except through the scalar product = k ■ q. Thus without loss of 
generality we can write f\ = fi{k, q, /i, r), and decompose /i in terms of a Legendre series, 

oo 
n=0 

*„(A;,g,T) = 2{-i)^ J-i^^ h{k, f^,q,r)Pn{fJ'), (23) 

with P„(/i) the Legendre polynomial of degree n. 

The advantage of this decomposition is manifest once we recognise that each moment ^„ 
can in fact be related to a bulk quantity of the neutrino fluid. From the zeroth moment we 
recover the energy density and pressure perturbations, 

while the velocity divergence and the anisotropic stress, 

d^q f q' 



{2ny 
■ d^q 
3"' J (27^3^ Ve; 



{Pu+Pu)a. = la-' f^,e(^-)\„ (25) 
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are linked to the first and second moments respectively. A set of evolution equations for the 
moments "ifn can now be constructed by integrating equation (20) as per equation (23): 



e d In g ' 



3e 3g 0? In g ' 

^n>2 = j^^^^^[n^n-l-{n+l)'^n+l]- (26) 

The is the so-called Boltzmann hierarchy. The multipole at which to truncate the hierarchy 
depends on the desired level of accuracy. In the Boltzmann code CAME [15], 1000 momen- 
tum bins and nmax = 6 appear to suffice for the computation of the matter power spectrum 
accurate to the percent level. 

3.4 Initial conditions 

Initial conditions for the perturbations are set by inflation (see, e.g., [14]). In the inflationary 
paradigm, a scalar field — the inflaton — ^whose dynamics is dominated by its potential energy 
drives an early phase of exponential expansion. During this expansion, quantum fluctuations 
on the scalar fleld are stretched to superhorizon scales {k ^ "H) and imprinted on the 
spacctime metric. At the end inflation, the inflaton releases its kinetic energy through particle 
production, a process known as reheating. When the kinetic energy of these particles come 
to dominate the energy budget, the universe enters into the radiation-domination era. 

Because the particles thus produced originate from the inflaton fleld, they must also 
inherit its superhorizon fluctuations. In the simplest scenario, the single-fleld inflation sce- 
nario, in which only one scalar fleld is responsible for both driving the exponential expansion 
and seeding the metric perturbations, the generic predictions for the initial perturbations 
are [11] 



1 2or^ / 2 \ 

a. = -(krf^^, *=(l + 5fl>. (27) 

where C is some constant determined by the amplitude of the fluctuations from inflation, and 

Ru = pv/ {pv + P7) during radiation domination. These are called adiabatic perturbations, 
and can be interpreted to mean that the number density ratios between the various particle 
species are the same everywhere. For the neutrino Boltzmann hierarchy, equation (27) implies 
that *o = -{K/Wk/dhiq), ^1 = -{ej^qk){dfo/d\Yiql and^^2 = -(a,/2)(rf/o/rfln g). 

Inflation does not predict exactly what initial values 5{k), ■?/'(^),and so on should adopt. 
Rather, it predicts that these quantities are random in nature, and their statistics are quan- 
tified by the n-point correlators. Of particular importance is the two-point correlator, e.g., 
(5(k)5(q)) = 5j:){k. + q)Ps{k), where Ps{k) is called the power spectrum of 6. Single-field 
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inflation models generically predict P{k) to be the only nontrivial statistics: all odd corre- 
lators are vanishingly small, while all even correlators can be constructed from P{k). These 
are known as gaussian initial conditions. Note that P{k) depends only on the magnitude 
and not the direction of k because of statistical homogeneity and isotropy. 

4 Neutrino mass and cosmological observables 

4.1 Observational probes and signatures of neutrino mass 

Massive neutrinos impact on cosmological observables such as the CMB and the LSS in 
two distinct ways. Firstly, neutrino masses in the sub-eV to eV range significantly alter 
the expansion history near the epoch of matter-radiation equality. Secondly, neutrino free- 
streaming affects the growth of structures at late times. In this section, I describe the 
manifestation of these effects in the CMB and the LSS as well as various observational 
probes currently used to detect these signatures of neutrino mass. 

4.1.1 Cosmic microwave background anisotropies 

The state-of-the-art instrument used to measure the CMB temperature and polarisation 
anisotropies is the Wilkinson Microwave Anisotropy Probe (WMAP), which measures the 
full-sky anisotropies with a resolution of 0.5 degrees [16]. Anisotropies on smaller scales have 
been measured by an array of ground-based and balloon experiments such as ACBAR [17], 
BICEP [18], QuAD [19], and ACT [20]. 

Because these fluctuations are mapped onto the surface of a sphere, it is convenient 
to decompose them in terms of spherical harmonics (cf. Fourier decomposition in flat 
space). For example, for the temperature fluctuations, 

r AT 
aJm = (-0'47r j dn Ylje, <t>)^{e, 0), (28) 

where 6 and 6 denote the angular position of a measurement, and dO, the solid angle. The 
statistics of the fluctuations are quantified at the lowest order by a two-point correlator 
{(if„fiji*n') = (^e('(^m,m'Cf^ , whcrc X.Y = T.E.B denotes temperature, and the E and B 
mode polarisation.^ Note that Cf^ has no dependence on the index m because of statistical 
homogeneity and isotropy. Figure 1 shows the temperature auto-correlation Cj-'^ for several 
different cosmologies. 

^Polarisation of the CMB is decomposed into a curl-free E-mode and a divergence-free B-mode. The 
latter is present at linear order only if there are primordial gravitational waves (tensor modes) or vector 
modes, and has yet to be detected. 
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Figure 1: The cosmic microwave background temperature anisotropy spectrum for a model 
with massless neutrinos and two models with massive neutrinos. Data are from WMAP after 
three years of observation [21]. 

The main signature of neutrino masses in the CMB comes about via the so-called early 
integrated Sachs- Wolfe (ISW) effect. Photon decoupling occurs at a temperature of T ~ 
0.26 eV. Shortly before this event, at T ~ 1 eV, is the epoch of matter-radiation equality. 
The transition from radiation to matter domination induces a nontrivial evolution for the 
metric perturbations and ip, and hence the photon geodesies. Because this transition occurs 
so close in time to the epoch of photon decoupling, remnants of this nontrivial evolution are 
preserved in the CMB anisotropies, especially in the height of the first acoustic peak of the 
CMB temperature power spectrum. This is the early ISW effect. 

If the absolute neutrino mass is roughly of order 1 eV, then the relic neutrinos also transit 
from being a fully relativistic particle species to a nonrelativistic one at T ~ 1 eV, thereby 
contributing to the early ISW effect. Figure 1 shows examples of the CMB temperature power 
spectrum for cosmological models with massive neutrinos, juxtaposed with the prediction of 
a model with massless neutrinos. This enhancement of the first peak enables us to put a 
constraint on the absolute neutrino mass scale from the CMB alone. 

4.1.2 Probes of large-scale structure 

The idea of "weighing neutrinos" with LSS probes hinges on two effects [22]. The first effect 
is the manifestation in the matter power spectrum of the background effect described above 
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Figure 2: The large-scale matter power spectrum for a model with massless neutrinos and 
two models with massive neutrinos. The data points are from the 2dF survey [24]. 

associated with the transition from radiation to matter domination. This effect is most 
important for sub-eV to eV mass. The second effect is free-streaming. Neutrinos are born 
with a large thermal velocity, which makes them difficult to capture into potential wells 
corresponding to wavenumbers larger than the so-called free-streaming wavenumber [23], 

even after they have become nonrelativistic. 
In other words, if we were to replace part of the CDM content with neutrinos, then free- 
streaming would cause the growth of the overall density perturbations at k > kf^ to be less 
efficient. 

In the large-scale matter power spectrum, these two effects combine to produce a sup- 
pression of power on small scales scaling roughly as AP/P ~ Sf^, where fu = Qu/^m is 
the neutrino fraction, according to linear perturbation theory. A detailed derivation can be 
found in reference [5] . Figure 2 shows this suppression for the case of J2 = 1 eV. Note 
that the suppression relative to the massless neutrino case depends primarily on the sum 
Yl,^v Exactly how the neutrino masses are distributed contributes only a minor effect. 

Several types of surveys of the large-scale matter distribution are available, each em- 
ploying different tracers of the matter density field. Currently, galaxy redshift surveys are 
the most powerful LSS probe, with the the Sloan Digital Sky Survey (SDSS) [25] being the 
largest to date. These surveys produce essentially galaxy catalogues, from which one can 
compute the galaxy number density at each point in space, and hence their n-point statistics. 
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The present generation of surveys observe at relatively low redshifts {z < 0.5) and at small 
length scales {k > 0.02 h Mpc""*^) just beyond the reach of CMB observations. 

Because galaxy surveys measure the galaxy number density fluctuations, in order to use 
these measurements to constrain cosmology, we must assume that the clustering of galaxies 
traces that of the underlying matter distribution up to a constant bias factor, i.e., 5gai = 
bSui{k). The value of b depends on the galaxy sample used, and is usually marginalised at the 
end of an analysis because it cannot be computed from first principle. This procedure works 
well on large scales, but is expected to fail on small scales. Indeed, analysis of the luminous 
red galaxy (LRG) sample of the SDSS already requires that we introduce some form of 
nonlinear modelling to account for the so-called scale-dependent bias at k > 0.1 h Mpc~^ [26- 
28]. However, because the nuisance parameters of the bias model are always marginalised 
at the end of the analysis, this procedure is effectively the same as discarding all (nonlinear) 
data beyond k > 0.1 h Mpc"^ [28,29]. 

Another LSS probe that has gained considerable attention recently is the Lyman-a (Lya) 
forest. In essence these are absorption features in intergalactic low density gases along the 
line of sight to a distant quasar. Since density fluctuations in these gases are expected to 
follow the underlying matter density fleld, measurements of the Lya flux power spectrum 
can be used to reconstruct the matter power spectrum on small scales at large redshifts 
{z > 1). However, unlike the case of the galaxy bias, the gas-matter bias is highly nontrivial, 
and shows strong dependence on the gas temperature and ionisation history, metal line 
contamination, and so on. So far various working groups have not reached a consensus 
regarding the imphcations of Lya for precision cosmology [30-33] . 

A number of other probes, such as weak gravitational lensing and cluster abundance, have 
also been explored to constrain cosmology. However, since these probes tend to rely strongly 
on nonlinear physics which is presently not very well understood, I defer their discussion to 
section 7. 

4.1.3 Other probes 

Standard candles and standard rulers are objects of known luminosity and known physical 
size respectively. Measurement of a standard candle's apparent magnitude allows us to 
deduce its luminosity distance. Similarly, the angular size of a standard ruler reveals its 
angular diameter distance. If the redshifts of a series of standard candles or rulers are 
known, then it is possible to establish the distance-redshift relation, and hence the low- 
redshift expansion history of the universe. The most well-known standard candles to date 
are the Type la supernovae (SN), while the position of the baryon acoustic oscillations (BAO) 
peak is a standard ruler frequently used in precision cosmology [34]. 

Another important measurement is the present expansion rate of the universe, quanti- 
fied by the present day Hubble parameter Hq. Its value — at least that in our immediate 
neighbourhood {z < 0.1) — has been estabhshed from observations with the Hubble Space 
Telescope (HST) to be Hq = 74.2 ± 3.6 km s"^ Mpc"^ [35]. 

Neutrino masses have little effect on the distance relations and the universal expansion 
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at low redshifts. Therefore these probes play no direct role in constraining the absolute 
neutrino mass scale. Nonetheless, they remain useful for removing parameter degeneracies. 

4.2 Cosmological model 

The simplest model required to account for all cosmological observations to date is the 
concordance fiat ACDM model. It makes the following assumptions. 

1. General relativity holds on all length scales. 

2. The large-scale spatial geometry of the universe is flat. 

3. The energy content of the universe comprises photons, whose energy density is fixed 
by the COBE FIRAS measurement of the CMB temperature and energy spectrum [7], 
3.046 families of thermalised massless neutrinos with temperature given by equa- 
tion (4), and unspecified amounts of baryons, CDM, and vacuum energy due to a 
cosmological constant. 

4. The initial conditions, i.e., the statistics and amplitude of the primordial scalar per- 
turbations to the FLRW metric, are set by the simplest, single-field infiation models. 
The perturbations are gaussian and adiabatic, and the presence of primordial vector or 
tensor modes is negligible. The power spectrum of density perturbations is described 
by a constant spectral index Ug and an amplitude As, i.e., Ps{k) — Agk"'", both of 
which are free to vary. 

In addition to these parameters, a minimal model of the CMB anisotropics requires an 
additional astrophysical parameter, the optical depth to reionisation r. In all there are six 
free physical parameters in the "vanilla" model of cosmology: 

9,bh'^,ilcBMh'^,^A,ns,As,r. (29) 

Spatial flatness implies fib + flcBU + = 1, from which we can deduce the reduced Hubble 
parameter h. 

Of course, these six parameters are not the only ones that can affect cosmology: there 
are many more. These other parameters are not included in "standard" analyses because 
they do not change significantly the goodness-of-fit, i.e., the best fit of the model does 
not improve by more than 3 or 4 with the inclusion of any one extra parameter. This 
improvement is deemed too modest by many practitioners to justify the inclusion of extra 
parameters, although it is unclear where one should draw the line. Indeed, some of these 
extra parameters, such as the neutrino mass, are physically extremely well motivated. 

The simplest neutrino mass limits can be obtained by including in the vanilla description 
a variable amount of neutrino dark matter characterized by the sum of the neutrino masses 
Y,^v- We call this the "vanilla -l-mj," model. 



17 



4.3 Present constraints and model dependence 

Using the effects discussed above, many authors have placed bounds on the sum of the 
neutrino masses. I discuss some of them below. 

4.3.1 Vanilla bounds 

Within the vanilla+mi/ framework, the 7-year data from WMAP (WMAP7) already limit the 
sum of the neutrino masses to Y^rrii, < 1.3 eV (95% credible interval) [13]. Inchiding small- 
scale CMB experiments in the analysis leads to virtually no improvement [36]. Incorporating 
distance information from BAO [34] tightens the bound to Y^m^ < 0.85 cV [36], while com- 
bining CMB with the SDSS LRG halo power spectrum (HPS) [27] gives i^u < 0.61 eV [36].^ 

Using also the HST measurement of Hq further strengthens the bound. For example, 
CMB+HST yields E^!. < 0.58 eV (95% C.I.) [36], while combining CMB+HPS+HST gives 
J2i^u < 0.44 cV [13, 36, 37]. Further improvement can be expected from the inclusion of SN 
data. Reference [38] finds < 0.28 eV (95% C.l.) for WMAP5+SN+BA0+MegaZ- 

-t-HST, where MegaZ is a SDSS photometric redshift catalogue and SN refers to the first 
year data set of SNLS [39] . However, one should exercise caution here regarding the use of 
SN data, for they appear at present to be plagued by unresolved systematic issues associated 
with the light-curve fitting methods [40] ; Discrepant light curve fitters have led to estimates 
of dark energy parameters that differ by as much as 2.5a. 

Several other LSS probes have also produced interesting constraints on J2^u, although 
they are arguably more strongly dependent on nonlinear physics than those quoted above. 
For example, reference [41] finds X^mj/ < 1.1 cV (95% C.I.) from weak lensing observations 
with the CFHTLS together with WMAP5. Another analysis combines the cluster mass 
function from weak lensing with WMAP3 to obtain J2 f^u < 1-43 eV (95% C.l) [42]. A recent 
variant uses the X-ray cluster abundance from the ROSAT All Sky Survey and the Chandra 
X-ray observatory with WMAP5+SN+BA0 to find Em^ < 0.33 eV (95% C.L) [43]. 

4.3.2 Extended models 

Cosmological models beyond vanilla-l-m,^ fall into two categories: those that "perturb" 
around vanilla-l-m,^, and those that represent a complete change of paradigm. 

Perturbations around vanilla Models that perturb around vanilla-l-m;^ include the in- 
troduction of (i) extended infiation physics, (ii) extra relativistic species, (iii) dynamical dark 
energy/modified gravity, (iv) non-flat spatial geometry. Each type of model is discussed be- 
low. 

(i) Popular extensions to the vanilla description of the primordial perturbations include 
a running spectral index (i.e., scale-dependent n^), the presence of primordial gravitational 
waves, and isocurvature modes (from, e.g., multi-field infiation). The latter two affect only 

^Despite their sensitivity to neutrino free-streaming, LSS measurements cannot constrain a 7-parameter 
model on their own, and must be used together with large-scale CMB observations. 



18 



the CMB anisotropies at low multipoles and are not directly degenerate with neutrino 
masses [44]. A running spectral index may mimic or offset to an extent the small-scale 
suppression in the matter power spectrum caused by neutrino free-streaming. However, any 
such fine-tuning necessarily leads to modifications of the large-scale CMB anisotropies and 
can therefore be constrained. Indeed, constraints on J2 from WMAP5 alone are already 
completely independent of these additional features [37] . 

(ii) Solutions to several outstanding particle physics puzzles predict the production of 
additional light particles in the early universe. Amongst these are light sterile neutrinos 
and hot dark matter axions. Such particles contribute relativistic energy density to the 
total energy budget, leading to iVeff > 3.046. As a CMB fit parameter, A^eff is highly 
degenerate with ifo, ^mh^ and indirectly X^mj,, so that in order to obtain constraints on 
X^mj, it is necessary to introduce additional data sets. The best available bound for the 
vanilla+m^+A^cff model is < 0.89 eV (95% C.I.), coming from CMB+HST+BAO [29]. 
The corresponding preferred value of N^q is 4.47lJ;74. 

(iii) Dynamical dark energy scenarios involve replacing the cosmological constant with 
a fluid whose equation of state parameter wde satisfies wde < —1/3 and may additionally 
be time-dependent [45]. More elaborate variants usually involve coupling between the dark 
energy and other matter components (e.g., [46]). Phenomenologically, Wde exhibits consid- 
erable degeneracy with the neutrino mass [47]. For a constant Wde, the best bound comes 
from WMAP7+HPS+HST, < 0.71 eV (95% C.I.) for vanilla+m^+WDE [13] (note that 
the CMB alone cannot constrain scenarios with a free wde)- This is to be compared with 
E?7ij/ < 0.44 eV for wde = — 1 using the same data sets. If in addition Aefr is allowed to 
vary, then the bound relaxes to J^n^-u < 1-16 eV (95% C.I.) [29]. 

Scenarios that modify general relativity at very large distances are primarily constructed 
to explain the observed late-time accelerated expansion of the universe in lieu of a cosmolog- 
ical constant. Phenomenologically they share many similarities with dynamical dark energy 
scenarios. 

An interesting issue concerns the possibility that certain coupled dark energy and mod- 
ified gravity scenarios can lead to scale-dependent clustering and thus mimic neutrino free- 
streaming [48]. This issue has yet to be explored in detail. 

(iv) Flat spatial geometry is one of the pillars of the inflationary paradigm; Relaxing the 
assumption of spatial flatness is perhaps the least theoretically well-motivated modifications 
to the vanilla model. Phenomenologically, however, a nonzero spatial curvature does appear 
to be degenerate with E^i/ [49], in the same way that w is degenerate with E^i^- 

Different paradigm Several cosmological models that differ radically from concordance 
ACDM have been proposed in the literature, including broken scale invariance [50] and 
inhomogeneous void models that seek to explain current cosmological observations without 
invoking a phase of late-time accelerated expansion. Although these models have found 
some success with subsets of the available data, it is generally difficult to reconcile the 
simplest variants with all data sets. The role of massive neutrinos has yet to be elucidated 
for inhomogeneous models. The model of [50], although incompatible with BAO and SN, 
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requires three species of neutrinos, each which has a mass of ~ 0.5 eV [51]. This value an 
be tested by the tritium /3-decay experiment KATRIN [52] . 

5 Nonlinear regime I: A/^-body simulations 

Linear perturbation theory has served the current generation of observational probes well. 
To advance further, however, future LSS probes will begin to observe on smaller length 
scales, where we expect the evolution of density perturbations to enter the mildly nonlinear 
regime. At the same time, scales that have already been observed will be probed with ever 
increasing precision. Thus in order to maximise the potential of future probes to constrain 
cosmology, we must go beyond linear order for accurate computations of the observables 
from theory. The length scales of interest here fall into the range k ~ 0.1 1 Mpc~^. 

To this end, two lines of approach, A^-body simulations and semi-analytic techniques, are 
currently being nvestigated in an attempt to bring the large-scale matter power spectrum 
under control at the percent level. The former is arguably the more definitive approach of 
the two and its validity generally extends to much smaller length scales than are possible 
with the latter. However, where they are valid, semi-analytic techniques do in general offer 
a significant reduction in computation time, and are in many ways more transparent. The 
two approaches therefore complement each other. 

In this section I describe the basic idea of how to compute nonlinear corrections to the LSS 
distribution, especially in the presence of massive neutrinos, by way of A^-body simulations. 
For a review of A^-body methods, see, e.g., reference [53]. Semi-analytic techniques will be 
discussed in section 6. Note that most current nonhnear investigations neglect the effects 
of baryonic physics, i.e., baryons are treated as a coUisionless matter component like CDM, 
while scattering processes responsible for the heating and cooling of gasses, star formation, 
and so on are ignored. The common conception is that these effects arc important only for 
k > 1 Mpc~^. Since we are concerned mainly with the mildly nonlinear scales /c ~ 0.1 ^ 
1 Mpc"\ it seems justified to ignore these effects. However, although baryonic physics is 
anticipated to be a minor player on the scales of interest, it is not inconceivable that some 
effect is present, even if only at the percent level [54]. Therefore, one should hear in mind 
that although the nonlinear techniques to be discussed below are all valid and correct in their 
own right, they probably represent only part of the full picture of nonlinear corrections. 

5.1 Basic idea of an A^-body simulation 

A^-body simulations are arguably the most definitive approach to nonlinear clustering to 
date. In this section, I described the basic building blocks of an A^-body simulation. 

5.1.1 Cosmological simulations 

Modern cosmological simulations are performed using a cubic "box" typically with a comov- 
ing side length of ~ 100 — )■ 1000 Mpc. This box is filled with a fiuid (or several fiuids) whose 
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phase space density fa evolves according to the Boltzmann equation (18) with the coUision 
term C[fa] set to zero. 

At z < 100, all particle species of interest are nonrelativistic. Rclativistic energy densities 
and the associated anisotropic stresses are negligibly small, so that the metric perturbations 
ip and (p are equal. Furthermore, we expect nonlinear evolution to be important only on 
subhorizon scales {k :$> %). These considerations justify the use of a pseudo- Newtonian 
description for the system, in the sense that apart from the Hubble expansion, all other 
aspects of gravitational physics are to be described by Newton's laws. In practice this means 

— = , — = -ama\/(p, (30) 

or rricea ar 

so that the coUisionless Boltzmann equation can now be written as 

^ + ^ ■ V/„ - am.V0 ■ ^ = 0, (31) 
OT mad aq 



where the Poisson equation 



a 



links the density perturbations to the gravitational potential, with VLa{T) = Pa (''")/Pcrit ('?")• 

The coUisionless Boltzmann equation expresses conservation of phase space density along 
characteristics {x(t), q(T)}. An fluid element drawn initially from {xj, qj} — > {xj + dx, q^ + 
d(\\ will move according to the characteristic equations (30). Thus tracking characteristics 
originating from all possible {xj, q^} is equivalent to solving the coUisionless Boltzmann 
equation. In practice one can only track a small sample of characteristics. This is then an 
A^-body simulation. 

The Boltzmann equation (31) imphes that fluid elements should be sampled from both 
coordinate and momentum space. However, since the emergence of the concordance ACDM 
model as the standard model of cosmology, a default cosmological simulation of CDM requires 
only that we sample from coordinate space. This is a a reasonable approximation; Already 
at linear order we have seen that the CDM fluid is treated as a single stream fluid described 
only by its density contrast and bulk velocity. The velocity dispersion is by deflnition zero. 

More precisely, we can take kinetic moments of the phase space density, 

/d?(l 
(27ra)^ '^"' ^^^^ moment) 

PaW = rria / — fa, (1st moment) 

Pacra,ij = rua [ ,^ , ^'^^^ fa " PaUiUj, (2nd moment) (33) 
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and of the Boltzmann equation (31), 



Sa + V ■[{! + 5a) Ua] = 0, (0th moment) 

Ua + Tivia + {via ' V)ua + V0 + —Vj{pa(Ta,ij) = 0. (ist moment) (34) 

Pa 

Simulating the CDM case simply consists of setting the stress term (Tc^ij to zero, thereby 
truncating the hierarchy at the first moment equation. The remaining zeroth and first mo- 
ments equations are none but the continuity and Euler equation in Eulerian fluid dynamics. 
Defining the convective derivative (i.e., the time derivative taken along a path moving with 
velocity Uq), D/Dt = d/dr + ■ V, we see that the Euler equation can be written as 
Dua/Dr + "HUq, = — V</). This equation is exactly equivalent to the characteristic equa- 
tions (30) if we identify = dx/dr. In other words, an element of CDM fluid moving with 
velocity u evolves according to Newton's second law of motion. 

Thus a recipe for simulating CDM emerges. Imagine a cubic box of comoving volume V 
filled with a fluid whose initial mass density is p . This box is divided into N blocks, with each 
block represented by a point particle of mass m = p/NV. These particles move imder one 
another's gravity according to equation (30), where Hubble expansion has been absorbed 
into the comoving coordinates. Periodic boundary conditions are imposed, meaning that 
one should envisage the box as one amongst an infinite lattice of identical boxes; Should a 
trajectory lead outside the box on one side, it would reenter from the opposite side. 

5.1.2 Initial conditions 

Simulations are usually started at high rcdshifts, typically 2; = 50 — )■ 100, where perturbation 
theory is still valid, and can be used to set the initial conditions of a simulation. The most 
popular implementation is by way of the Zel'dovich approximation and its variants. 

Imagine a simulation box spanned by a 3D grid, where we place on each grid point point 
particles of identical masses. The particles thus positioned represent a uniform fluid. To 
factor in the initial density contrasts and bulk velocities, we assign to each particle a random 
initial displacement s(r) from the grid point r and an initial velocity vector v as dictated 
by the (perturbative) clustering statistics at the time the simulation is initiated r^. In the 
Zel'dovich approximation, the displacement vector s is given by V • s = —S{r,Ti)/D^{Ti), 
where -D+(Ti) is the linear growth function at the initial time, 5(r, r,) is a random number 
drawn from a gaussian distribution with variance determined by the linear power spectrum 
P{k, Ti), and V x s = 0. Thus, the initial position of the simulation particle and its velocity 
are now given by 

x-r + L>+(T,)s(r), v = — s(r). (35) 

or 

The randomly displaced particles now represent an inhomogeneous fluid whose clustering 
statistics respects the predictions of linear perturbation theory; 

The Zel'dovich approximation is technically the first order solution to the fluid equa- 
tions (34) in Lagrangian perturbation theory. In modern A^-body simulations, the displace- 
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ment vectors and initial velocities are generally computed from second order Lagrangian 
perturbation theory. 

5.1.3 Gravity solver 

The heart of an A^-body code is the gravity solver, which computes the gravitational force 
felt by each simulation particle. Direct summation of the forces between each pair of particles 
may be intuitive. However, because the number of operations scales as N"^ where N is the 
number of particles, this approach quickly gets out of hand in cosmological simulations in 
which N often exceeds a million. 

Two broad classes of gravity solvers are employed in modern simulations: the Particle- 
Mesh (PM) code, and the Tree code. The PM code is particularly important for cosmological 
simulations because of the ease with which one can implement periodic boundary conditions. 

Peirticle-Mesh code The PM code uses a grid or mesh on which to evaluate the grav- 
itational forces. The presence of a mesh enables the use of discrete Fourier transforms via 
the Fast Fourier Transform algorithm, a pivotal routine in this class of gravity solver. Here, 
for the purpose of evaluating 0(x), the particles are envisaged as clouds and then smoothed 
on to the grid points. The total mass density p(x) on each grid point is then a sum of each 
particle's contribution. This density is then Fourier transformed into p(k), which allows for 
an algebraic evaluation of the gravitational force (in Fourier space) J^[V0] via the Poisson 
equation (32), i.e., J-'[V0] = ik0 = —ihAirGa'^p/k'^. This value is inverse Fourier transformed 
back to real space to give V<6 on each grid point, and then interpolated to the particle posi- 
tions to be fed into the characteristic equations (30). The total number of operations scales 
as Nm log Njn, where Nm is the number of grid points, significantly less than that required in 
a direct summation scheme. The disadvantage, however, is that the force resolution becomes 
unrehable on length scales close to and below the grid spacing. 

Tree code To illustrate the principle of the Tree code, take as an example three particles 
A, B and C, where B and C are separated by a distance s, and d denotes the distance from A 
to the centre-of-mass (c.o.m.) of B and C. To evaluate the combined force of B and C on A, 
the Tree algorithm stipulates that we treat B and C as one particle centred on their c.o.m. 
if their angular separation satisfies s/d < 6, where 9 is some small number of order unity. If 
s/d > 6, then B and C must be treated as two separate particles. This idea of lumping two 
particles together according to their angular separation can be easily generalised to the multi- 
particle case, and has the advantage that the total number of operations scales as N log N. 
A popular implementation is the Barnes-Hut algorithm [55]. However, periodic boundary 
conditions are not intrinsic to the Tree algorithm as they are to the PM code. Nonetheless, 
the Tree algorithm is often used in combination with the PM code in cosmological simulations 
to overcome the small-scale shortcomings of the latter, such as in the publicly available A'^- 
body code GADGET-II [56]. 
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5.2 A/^-body simulations with massive neutrinos 

Simulation with massive neutrinos have many parallels with pure CDM simulations, but 
are complicated by the fact that neutrinos always come with large thermal velocities. In 
this section, I describe how large neutrino thermal velocities translate into complications for 
A^-body simulations and review some recent efforts towards meeting these challenges. 

5.2.1 Thermal neutrino motion 

Simulating massive neutrinos is nontrivial because of their large thermal velocities. Recall 
that the unperturbed neutrino phase space density is described by a relativistic Fermi-Dirac 
distribution. For a 1 eV neutrino, equation (6) shows that the average thermal velocity is 
easily 15% the speed of hght at z = 50 when we start a simulation, even though the bulk 
velocity of the neutrino fluid is much less than the Hubble flow. 

To simulate (nonrelativistic) massive neutrinos, we employ the same equations of motion 
as those for simulating a CDM fluid. Relativistic corrections are generally too small for 
the desired level of accuracy to be of concern. The only difference between neutrinos and 
CDM arises when the initial conditions are implemented: in addition to sampling from 
coordinate space, it is necessary to sample from momentum space for a neutrino fluid to 
account for its velocity dispersion. A simple implementation introduced in reference [57] uses 
the Zel'dovich approximation to generate an initial random bulk velocity for the neutrino 
fluid at each grid point. Added to this bulk velocity is a thermal component consisting of a 
random thermal speed drawn from the relativistic Fermi-Dirac distribution and a random 
unit vector specifying the direction of the thermal contribution. 

Figure 3 shows the large-scale matter power spectrum for several massive neutrino cos- 
mologies from the simulations of reference [57], including the aforementioned neutrino ther- 
mal velocity dispersion (see also reference [58]). More precisely, this figure shows the sup- 
pression of power induced in the matter power spectrum when we replace part of the CDM 
content with massive neutrinos. Importantly, although linear perturbation theory predicts 
a relative suppression of AP/P ~ 8fu in the matter power spectrum on small scales due to 
neutrino free-streaming, the A^-body results in Figure 3 indicate that the suppression in fact 
scales more like AP/P ~ 10/,^! 

5.2.2 Approximate schemes 

Because of the additional random sampling involved, simulations including massive neutri- 
nos can be plagued by an increased level of shot noise unless the sampling is adequate. Also, 
the neutrinos' large thermal velocities mean that the time steps taken by the simulation 

code must be appropriately small in order to accurately track the trajectories. Starting 
the simulation at a lower rcdshift when the neutrinos have slowed down somewhat allevi- 
ates these problems to a good extent but at the expense of reduced accuracy.^ All in all, 

high initial redshift is required if one is to compute the absolute matter power spectrum to percent 
accuracy. However, if one is interested merely in the relative suppression due to neutrino free-streaming, 
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Figure 3: The large-scale matter power spectrum for massive neutrino cosmologies relative to 
the case with massless neutrinos from A^-body simulations. The solid/green lines represent 
the predictions from linear perturbation theory, while the colourful dotted/dashed lines 
are the results from full-scale nonlinear simulations for different neutrino masses {J2 ^u) 
indicated on the plot. Figure reproduced from reference [57]. 



simulating massive neutrinos can be very computationally intensive and time-consuming, so 
approximate schemes are worth consideration. 

Grid-based linear neutrino perturbations As the name implies, in this scheme the 
neutrino perturbations are assumed to remain at linear order and evolved on each grid point 
in a PM code using a linear Boltzmann code such as CAMB [15], while the CDM fluid is given 
a particle realisation. The neutrinos' contribution to the total mass density is then added 
to the CDM's contribution on each grid point, so that the gravitational force felt by each 
CDM particle can be evaluated via the usual PM machinery. This scheme was explored in 
reference [59] and found to be accurate to the percent level for J2 < 1 eV compared with 
a full-fledged simulation. However, the reduction in the simulation run time can be as much 
as a factor of ten. Note that in reference [59] nonlinear corrections to the CDM perturbations 

then starting a simulation as late as z = 4 does not affect the results [57]. 
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are not fed back into the evolution equations of the hnear neutrino perturbations, which may 
be a point to improve upon. 

Hybrid scheme A more adventurous hybrid scheme was introduced in reference [60]. 
Here, the neutrino perturbations for aU momenta are initiaUy tracked with hnear perturba- 
tion theory, a reasonable approximation given that at high redshifts the neutrino thermal 
speeds are so large that neutrino clustering must be minimal. As the neutrinos slow down at 
low redshifts. beginning with the low momentum states, those neutrinos with a momentum- 
to-temperature ratio below some predefined threshold are converted into A^-body particles 
upon crossing.^ 

This scheme avoids the initial shot noise problem, but at the same time allows us to 
track the neutrino nonlinearities fully at late times down to very small length scales. It is 
thus particularly well-suited to the study of such small scale problems such as halo density 
profiles and the halo mass function. The latter counts the number density of virialised 
structure ("halos") per unit mass at a given redshift, and was computed for massive neutrino 
cosmologies for the first time via A^-body simulations in reference [61]. Figure 4 shows the 
halo mass function for several cosmologies. 

6 Nonlinear regime II: Semi-analytic methods 

Numerical simulations are arguably the most definitive way to date to treat the nonlinear 
clustering problem. However, time, they can also be very computationally demanding and 
time consuming, especially for cosmologies beyond the ACDM model. For this reason, a 
number of semi-analytic schemes have been investigated for the calculation of nonlinear 
corrections to the matter power spectrum and other LSS observables. I describe two below. 

6.1 Higher Order Perturbation Theory 

As the name implies, the higher order perturbation theory approach consists of going beyond 
linear when solving the equations that describe the evolution of the density perturbations. 

6.1.1 The equations 

The starting equations in this approach are the fluid equations (34) and the Poisson equa- 
tion (32) already introduced earlier (for a detailed review of perturbative methods, see, e.g., 
[62]). Because the main constituent of the matter content is CDM, it makes sense to con- 
sider first nonlinearities in the CDM fiuid, while keeping the neutrino perturbations at linear 
order. This is the approach taken in references [63,64], and is a reasonable approximation 

^In practice, owing to the Eulerian nature of Boltzmann codes like CAMB and the Lagrangian nature 
of iV-body simulations, a bin by bin conversion is a nontrivial accounting problem. A simple solution that 
works well is to convert a large number of momentum bins at the same time. See discussion in reference [60] . 
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Figure 4: Halo mass function for cosmologies with the indicated neutrino masses [J^^i^)- 
The 0.0 eV case corresponds to the ACDM model. Figure reproduced from reference [61]. 

supported by A^-body simulation results [59] for neutrino masses not exceeding J2 '^u ~ 1 eV. 
There have also been attempts to model the nonlinearities in the neutrino fluid [65]. 

Cold dark matter is pressureless at linear order so that a^^Jj = 0. At very nonlinear scales, 
however, we do expect pressure and anisotropic stress to be generated. However, because we 
are interested only in the mildly nonlinear regime (fc ~ 0.1 — )■ 1 Mpc~^), these factors are 
probably not very important. Therefore the first approximation is to assume that (Tc,ij = 
at all orders. It follows that in the absence of a finite Cc^ij, no vorticity = V x can be 
generated for the CDM fluid. 

With these assumptions, we can now rewrite the continuity and the Euler equations in 
terms of the only remaining degrees of freedom, the density contrast 6c and the velocity 
divergence = V • u^. Since we are working with perturbation theory, it is convenient to 
rewrite these equations in Fourier space, 

6c(k, r) + 6'c(k, r) = - J d^qi(fq25D(k - qi2)«(qi, q2)^c(qi, r)5c{q2, r), 
^^(k,r)+m(k,r)-A;20(k,r) = 

J d^qid^q2SD{k - qi2);5(qi, q2)^c(qi, r)^,(q2, r), (36) 

where qi2 = qi + q2, and the vertices a(qi,q2) = qi2 ■ qi/g? and /3(qi,q2) = Qui^i ' 
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qa) / (2^1 Q'l) encode the coupling between Fourier modes. 

At this point, we define a doublet ip = [5c, —Oc/T^Y^^ ^-nd a new time variable s = In a. 
Then the equations of motion (36) can be recast into a compact form, 

ds^aiK s) + n„,,(k, s)v36(k, s) = j d^qid^q2'yabciK Qi, q2)¥'6(qi, s)<^c(q2, s), (37) 

where 



U{k,s) 



-1 



;[^^c(^) + Ea^c^^a(^)ff^] 1+^'^^ 



(3J 



and the vertex function ^yabc is zero except for 

7i2i(k,qi,q2) = 7ii2(k, qa, qi) = ^^(k - qi2)^a(qi, qa), 

7222(k,qi,q2) = 5i3(k - qi2)^(qi, q2). (39) 

Let us call equation (37) the "master" equation, as we will return to it below. 

In standard perturbation theory, we solve equation (37) using a perturbative series, 
</?(k, s) — J2'^=i </7^"'^(k, s). The equation of motion for the n-th order perturbation is therefore 

a,(^i")(k,s) + n„,(k,5)^i")(k,s) = 

/ d%d%^,U^, qi, q2) X: 'Pt'"\^u s)^'r\^, s), (40) 

m=l 

which has the formal solution 

ip^:\k,s)=g,,{Ks,s,)^t'\Ks.) (41) 

+ [d\ [d% fds'ga,{K S, s'hSJ: V^i^-^Hqi, ^>'r\^2, S'). 
J J Jsi 

Following field theory language, 5'ab(k, s, s') is known as the linear propagator. For a matter- 
dominated universe, i.e., Oc(s) = 1) 5*06 = 5'a6(S) s') is independent of A;, and takes the form 

^') = -^( 3 2] ^[ 3 -aj- (^2) 

The first term is a growing solution; The astute reader will have recognised that this is but 
the linear growth function for a matter-dominated universe. The second term is a decaying 
solution. 

The formal solution (41) in its present form has little practical use. To make it more 
user-friendly, we introduce some simplifications: (i) We assume that initially all higher order 
terms are vanishingly small; only (^'^^-'(k, Sj) is nonzero; (ii) Linear solutions are always in 
the growing mode. In other words, whenever we encounter a linear term </7^^^(k, s), we make 
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the approximation </'^^^(k, s) — **</?^^)(k, Sj); and (iii) The initial time of the time integral 
is in the infinite past, i.e., = — oo. Then, an order by order evaluation shows that 



(^i")(k, s) = J d\ ■ ■ ■ d%Sn{k - qi...„)gi")qi, . . . , c^r^)ip^\c^^, s) ■ ■ ■ (^S^^Qn, s), 
where the kernels are given by 

n-l 

Qi"^qi, ■ ■ ■ , Qn) = C^l? 7bcd(k, qi...m, Qm+l-n) 

m=l 

xg(™)(qi,...,qjg('^-™)(q^+i,...,q„), 

with 



An) _ 
^ab — 



1 



(2n + 3)(n- 1) 



2n + 1 2 
3 2n 



Note that in older literature, the kernels are denoted as F„ = Q^"^ and Gn = ^2"^- 



(43) 



(44) 
(45) 



6.1.2 The power spectrum 

We arc interested in the clustering statistics of matter as quantified by the matter power spec- 
trum ((5m(k, r)5m(q, t)) = ^^(k + q)Pm(^5 t), where the matter density contrast ^^(k, r) = 
(1— /,^)(5c(k, r)+/i,5(k, r) counts both contributions from CDM and neutrinos. Thus Pm{k, r) 
is a linear sum of the CDM and the neutrino auto-correlation spectra Pc{k, r) and Pu{k, r), 
and their cross-correlation spectrum Pcu{k,T), 

Pm{k, r) = (1 - U)'Pc{k, t) + 2(1 - U)UPUk, t) + f,P,{k, r). (46) 

Because we are keeping the neutrino perturbations at linear order, we need only to consider 
nonlinear corrections to Pc{k,T) and Pcu{k,T). 

As with the density perturbations, the power spectrum can also be expressed as an infinite 
sum. For the CDM auto-correlation, for example, we find 

P,{k, t) = Pl\K r) + P^iK r) + Pl\K r) • • • , (47) 

with 

5r,{\^ + c^)Pl\k,T) = (5«(k,r)5«(q,r)), 
5D(k + q)Pf(A;,r) = (5f (k, r)5p)(q, r)), 

5n{\^ + q}Pl\KT) = (5W(k,r)5f (q,T)) + (5f (k,r)5W(q,r)). (48) 

Here, the leading order term PJ^(A;,t) is the usual hnear power spectrum, while P'^'^ -\- P^^ 
are collectively known as the one-loop correction. Similarly for the CDM-neutrino cross- 
correlation, 

P,,{k, r) = P^'ik, t) + P^^ik, r), (49) 
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with 5D(k + ci)P^^{k,T) = (5(=^)(k,r)5(^)(q,r)). The "22" term is absent here because the 
neutrino perturbations are assumed to be hnear. 

Given the solutions (43), we can readily evaluate the higher order terms: 



where the subscript "(s)" denotes symmetrised (under the exchange of the wave-vectors 
in the argument) versions of the kernels given in equation (44). At this point the astute 
reader may question how a set of kernels originally derived for a pure CDM cosmology can 
be applied in this case of mixed CDM+neutrino cosmology, especially given that the hnear 
propagator in this case is /c-dependcnt. Corrections to these expressions have in fact been 
derived for this very case of mixed CDM+ncTitrino cosmology [64]. However, for the same 
small neutrino masses required to validate the linear neutrino perturbation assumption, the 
corrections at A; ~ 0.1 — )■ 1 Mpc~^ are very small — generally no more than a percent. In 
the same vein, the one-loop corrections (50) are often also extended to ACDM cosmologies 
without further modification [62]. 

Figure 5 shows the one-loop corrected matter power spectrum for several massive neutrino 
cosmologies relative to the ACDM case. The main feature is an enhanced suppression on 
small scales due to neutrino free-streaming: the suppression is stronger than is predicted 
from linear perturbation theory, which is qualitatively consistent with the outcome of N- 
body simulations in figure 3. 

6.2 Resummation/Renormalisation Group 

The time-flow renormalisation group (RG) method is a simple extension of the perturbative 
method introduced above, but incorporates resummation of certain classes of perturbative 
terms to all orders [66]. Aside from its simplicity and ease-of-use, the time-flow RG method 
has the advantage that it has been formulated for a broad class of cosmological models, 
including those with massive neutrinos. Other resummation approaches — of which there are 
many — are in their present formulations strictly speaking limited to ACDM or pure CDM 
type models. 

Recall the master equation (37) for Lp. Ultimately, we are not interested in the density 
contrast or velocity divergence per se, but rather the n-point statistics such as the power 
spectrum {ipa(k, s)(pb{q_, s)) = 6d0<- + (i)Pab{k;s). Therefore, it makes sense to construct 
directly equations of motion for the n-point correlators from the master equation (37). For 
the two-point correlator, we find 




(50) 




(51) 
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Figure 5: Relative differences between the total matter power spectra for a pure ACDM 
cosmology and three models with massive neutrinos, with = Q^/Qm = 0.1 (red/solid), 
0.05 (blue/dotted), and 0.01 (green/dash) at 2; = 1 (left) and z = 3 (right). Thick lines 
indicate results including the one-loop correction, while the linear results are represented by 
the thin lines. The three vertical lines indicate the maximum k values at which the linear 
and the one-loop corrected matter power spectra are accurate to better than 1% and 5%. 



Nonlinear coupling leads to a three-point correlator {ipafcy^d) to appear in the expression, 
which requires that we write out the corresponding equation of motion. However, nonlinear 
coupling necessarily introduces a four-point correlator into the equation of motion for the 
three-point correlator. Thus a pattern emerges: the presence nonlinear coupling gives rise to 
an infinite hierarchy of differential equations for the n-point correlators, and this hierarchy 
must be truncated if we wish to produce a result. 

In the time flow RG scheme, the hierarchy is truncated by simplifying the four-point 
correlator. In general the four-point function is given by 

where we have implicitly assumed (ipa) = (i.e., the mean density perturbation and mean 
velocity divergence are zero). The connected piece, denoted by the subscript "c", can be 
written as {(pa(pb'-PcVd)c = 5D(k + q + p + r)Q„fecd(k, q, P, r, s), where Qatcd is the trispectrum. 
We set Qabcd to zero. 

Thus, expressed in terms of the power spectrum Pab{k; s) and the bispectrum Babc{^, q, r; s), 
where {(faO^, s)(pb{(i, s)v?c(i", s)) = ^^(k -|- q -|- r)Babc{^, q, r; s), the final set of equations to 
be solved is [66], 

dsPabO^; s) = -n„c(k; s)Pbc{'k; s) - Ubcik; s)Pac{'k; s) 

+ J d^q [7acd(k, -q, q - k)BbcdiK -Q, Q - k; s) 

-h76cd(k, -q, q - k)Bacd{K -q, q - k; s)]. 
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dsBabciK -q, q - k; s) = -n„<i(k; T)Bdbc{k, -q, q - k; s) 

-n6<i(k; s)BadcO^, -q, q - k; r) - Ucd{k; r)Babd{K -q, q - k; s) 
+2[7ade(k, -q, q - k)Pbd{<i; s)-Pce(k - q; s) 

+76cie(-q, q - k, k)Pc(i(k - q; T)P„e(k; s) 

+7c<ie(q - k, k, -q)P„d(k; r)Pbe(q; s)]. (53) 

Reference [67] applies these equations to the massive neutrino problem, again assuming that 
the neutrino perturbations remain at linear order, and demonstrate that they reproduce the 
absolute matter power spectra from the A^-body simulations of reference [59] to a few percent 
accuracy up to A; ~ 1 ^ Mpc~^ at high redshifts {z > 2). This is superior to the one-loop 
perturbativc result, which has the same accuracy only up to A; ~ 0.4 h Mpc~^ at the same 
redshifts. Interestingly, however, the time-flow RG and the one-loop pcrturbative approaches 
perform equally well when it comes to predicting the relative suppression due to neutrino 
free-streaming at 2; > 2. Here, percent accuracy can be claimed up to /c ~ 1 /i Mpc~^. See 
figures 2 and 3 of reference [67] . 

7 Future observations 

7.1 Cosmic microwave background 

The Planck spacecraft, launched in May 2009, is expected to provide measurements of 
the CMB temperature fluctuations limited only by cosmic variance up to a multipole of 
i ~ 1000 [68]. Cosmic variance-limited observations of the E-mode polarisation are ex- 
pected to be possible up to i of several hundred. For the neutrino mass measurement, these 
improvements imply that Planck alone will have a sensitivity to J2 as good as 0.2 eV for 
the vanilla-l-m,^ model (e.g., [69]).^ Even for more complicated model frameworks including 
a variable and dynamical dark energy, we can anticipate (yij^mv) ~ 0.4 from Planck 
alone [69]. Importantly, notwithstanding the improved sensitivities, linear perturbation the- 
ory will remain applicable to CMB physics — no other observation will take its place as the 
most robust probe of cosmology in the near future. 

7.2 Large-scale structure 

Galaxy redshift surveys over the next decade will probe the LSS at higher redshifts than 
before (up to 2; ~ 4). The advantage is two- fold: firstly, nonlinear effects are less pronounced 
at early times, meaning that it will be possible to observe up to A; ~ 1 Mpc~^ with only 
mildly nonlinear corrections (however, scale-dependent galaxy bias will still be an issue [70]). 
Secondly, going to higher redshifts will enable us to observe a larger volume V", thereby 
reducing the sampling uncertainty which scales as AP/P ~ {Vk^)~^^'^. The HETDEX, 

^AU forecasted sensitivities given in this section refer to the vanilla-|-mi^ model unless otherwise indicated. 
Sensitivity is defined as the 68% upper limit an observation can place on J2 if ti"ue X) is zero. 
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WFMOS, and BOSS surveys are envisaged to push a{Yl rrii) down to 0.1 eV when combined 
with Planck, while with JDEM or EUCLID, a(Em^) ~ 0.05 eV may be possible [71-76]. 

Concerning Ly-a, barring uncertainties in the gas-bias relation, the Keck and/or VLT 
spectra of the Ly-a absorption lines have the potential to yield ai^my) ~ 0.05 [77] when 
combined with Planck. 

Perhaps more interestingly, new kinds of observations utilising different techniques to 
probe the large-scale structure distribution will also become available. I describe some 
examples below. 

7.2.1 Weak gravitational lensing of galaxies 

Light rays from distant galaxies are bent by matter density perturbations between the source 
galaxies and the observer, thereby inducing distortions in the observed images of the source 
galaxies. By measuring the angular correlation of these distortions, one can probe the 
clustering statistics of the intervening matter density field [78] . 

The current generation of weak lensing surveys are already producing interesting con- 
straints on the neutrino mass (see section 4). Future dedicated lensing surveys will probe 
higher redshifts (up to 2; ~ 3) with full sky coverage. Furthermore, all these surveys will 
provide photometric redshift information on the source galaxies, allowing for the binning of 
galaxy images by redshift and hence tomographic studies of the evolution of the interven- 
ing large scale structure. The Large Synoptic Survey Telescope (LSST) [79] combined with 
Planck offers a sensitivity of (yiJZi^v) ~ 0.04 eV using tomography. Similar sensitivity is 
expected for DUNE/Euclid [80-85]. 

Dominant systematics include uncertainties in the photometric redshift measurements, 
which are crucial for tomographic studies. Image distortions arising from the atmosphere 
and imperfections in the instrumentation must also be accurately known. Because future 
lensing surveys will derive most of their constraining power at nominally nonlinear scales 
k > 0.1 Mpc~^, uncertainties in our predictions of the nonlinear power spectrum must be 
controlled to the percent level. Baryon physics will also be important here. The study of [54] 
finds that baryon physics can contribute an uncertainty of up to 10% at multipole £ > 1000 
corresponding to physical scales k of less than a factor of ten beyond the linear regime. 

7.2.2 Weak lensing of the CMB 

Light from the last scattering surface can also be lensed by the foreground matter [86]. 
Indeed, with Planck one can hope to probe J^f^u with sensitivity of 0.1 eV via CMB lensing 
[69,87,88]. The main advantage of CMB lensing over galaxy lensing is that because the 
source lies at such a high redshift {z ~ 1000), the lenses probed are mainly those lying 
at 2; ~ 4, where the density perturbations are still evolving at linear order. The drawback, 
however, is that tomographic studies are not possible because there is only one last scattering 
surface. Another potential problem is that extraction of the lensing signal hinges on the 
non-gaussian nature of the signal. Other sources of non-gaussianity, whether primordial or 
secondary, may be an obstacle. 
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7.2.3 Cluster mass function 

As shown in figure 4, the halo mass function is sensitive to the absolute neutrino mass scale, 
because the formation of virialised objects is directly linked to the amplitude of the initial 
density perturbations. The halo mass function is also known as the cluster mass function, 
since it is the observation of galaxy clusters that ultimately allows us to identify virialised 
objects in the low redshift universe. 

As discussed in section 4, current X-ray observations of several hundred clusters have 
already produced interesting limits on the neutrinos mass. Future optical surveys such as the 
LSST are expected to yield cluster numbers in the thousands, thereby pushing to sensitivity 
down to cr{Y^m^) = 0.04 eV when their data are combined with data from Planck [89]. The 
main systematic challenge to these observations is how to determine the mass of the observed 
object. For X-ray observations, for example, the X-ray temperature is a proxy for the cluster 
mass; the relation between the two must be calibrated either through numerical simulations 
or through other observational means such as gravitational lensing effects. 

A variant of these methods uses the Sunyaev-Zel'dovich effect to probe the cluster mass 
function [90]. However, the expected a{J2'm'u) is only ~ 0.28 eV. 

7.2.4 21 cm 

Neutral hydrogen atoms can be found in two energy states depending on the alignment of 
the proton and the electron's spins. A spin flip may occur taking the atom from the high 
to the low energy state by emission of a photon with a wavelength of 21cm. Conversely, 
absorption of a 21cm photon brings the atom from the low to the high energy state. Neutral 
hydrogen is the main constituent of the baryonic content of the universe and its distribution 
in space at high redshifts is expected to trace the underlying dark matter density field. Thus, 
a measurement of the spatial fluctuations in the 21cm brightness temperature could in 
principle be a powerful way to map the LSS distribution [91]. 

Current and upcoming radio arrays such as LOFAR, MWA, and SKA are designed pri- 
marily to study the physics of reionisation at redshifts 6 < z < 12 when neutral hydrogen 
reside in the intergalactic medium. This is less than ideal for the purpose of mapping the 
LSS distribution, since fluctuations in in this epoch depend not only on the matter density 
perturbations, but also on the process of reionisation, which is probably a spatially inhomoge- 
neous process. If reionisation physics can be accurately modelled, then it should be possible 
to reach (rlJ^m^) ~ 0.02 eV using a combination of Planck and SKA [92, 93]. Gravitational 
lensing of the 21 cm signal could also yield similar sensitivities [94,95]. Innovative experi- 
mental design such as the proposed Fast Fourier Transform Telescope (FFTT) [96] could even 
yield 0.003 eV, which would enable one to probe even the individual neutrino mass [92, 93]. 
However, if no accurate modelling of reionisation is available, then the sensitivities to J2 f^v 
would significantly degrade, yielding aij^'^v) ~ 0.02 eV from Planck-|-FFTT [92]. 

Note that, with present technology, it is not possible to observe at higher redshifts to 
skirt the issue of reionisation. From the perspective of cosmological parameter estimation, 
observations at lower redshifts when neutral hydrogen forms dense clumps in galaxies offer 
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no real advantage over conventional galaxy redshift surveys. 

8 Conclusions 

Precision cosmology offers an interesting way to probe many aspects of fundamental physics. 
Neutrino properties are an excellent example. By using probes of the large-scale structure 
and of the cosmic microwave background anisotropics, we can already deduce that the sum 
of the neutrino masses must be less than about 1 eV regardless of the precise details of the 
cosmological model. This limit is already better than those from laboratory experiments. 
Future cosmological observations will perform even better: in the most optimistic case, a 
sensitivity of clJ^u) ~ 0-04 eV may be achievable with a combination of Planck, galaxy 
redshift surveys, weak gravitational lensing observations, and so on. 

Prom a theoretical perspective, the key issue is to maintain control over the nonlinear 
physics: specifically nonlinear evolution of the density perturbations and the galaxy bias. 
To this end, there have been recent developments in the computation of the nonlinear power 
spectrum and other observables to percent level accuracy by numerical and semi-analytical 
means. Understanding the nonlinear issue will aUow us to use future observations to their 
maximum potential, and hopefully to realise the goal of measuring the absolute neutrino 
mass scale with precision cosmology one day. 
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